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I.  OBJECTIVES 

Project  objectives  are  the  development  of  more  optimal  mechanics  approaches  for  textile 
composite  design  and  failure  analysis.  Tailoring  of  the  textile  composite  microstructure  is  one  of 
the  most  pressing  research  issues  in  textile  composite  design.  The  textile  pre-forming  process 
determines  the  microstructure  of  the  textile  preform.  Preform  microstructure  determines  textile 
composite  micro-stress  distribution.  Development  of  a  numerical  approach  that  facilitates 
establishment  of  relations  between  textile  microstructures  and  textile  processes  is,  therefore, 
critical.  In  this  project,  two  new  numerical  methods  are  developed.  The  first  is  a  digital  element 
simulation  approach  for  textile  mechanics.  It  enables  simulation  of  the  textile  process  as  well  as 
simulation  of  textile  preform  deformation.  As  a  result,  detailed  knowledge  of  the  textile  preform 
microstructure  becomes  obtainable.  The  second  method  developed  in  this  project  is  a 
heterogeneous  element  method  for  the  micro-stress  analysis  of  textile  composites.  Formulation 
of  a  heterogeneous  element  guarantees  that  both  the  equilibrium  conditions  and  the  continuity  of 
displacement  at  the  interface  are  satisfied.  It  also  allows  for  interface  stress  and  strain  jump. 
Because  the  formulation  reflects  the  actual  physical  situation  at  the  interface,  it  provides  a  much 
more  accurate  result  than  conventional  approaches  given  the  use  of  the  same  mesh.  This  project 
serves  the  interest  of  the  Air  Force  Research  Laboratory  because  it  contributes  to  the  goal  of  its 
on-going  project  “Micro-mechanical  Failure  Criteria.”  The  contact  persons  at  the  Air  Force 
Research  Laboratory  are  Senior  Scientists,  Dr.  Nick  Pagano  and  Dr.  Ajit  Roy. 

II.  DIGITAL  ELEMENT  ANALYSIS  IN  TEXTILE  MECHANICS 
2.1  Basic  Concept  of  Digital  Element  Analysis 

Digital  element  analysis  is  a  new  numerical  tool  developed  in  this  project  for  textile  mechanics. 
It  can  be  used  for  textile  process  design  and  fabric  deformation,  strength  and  failure  analysis. 

In  digital  element  analysis,  each  fiber  or  yam  is  modeled  as  a  frictionless  pin-connected  rod 
element  chain.  These  rod  elements  are  defined  as  “digital  elements.”  Contacts  between  fibers  or 
yams  are  modeled  by  contact  elements.  Textile  processes  and  fabric  deformation  are  formulated 
as  a  non-continuum  mechanics  problem  with  boundary  conditions.  A  procedure,  similar  to  finite 
element  analysis,  is  adopted  to  derive  yam  movement  during  textile  processes  or  fabric 
deformation. 

Three  critical  concepts  suffuse  digital  element  analysis:  digital  chain,  contact  between  digital 
chains,  and  yam  assembly. 


2 


2.1.1  Digital  chain 

A  digital  chain  is  composed  of  many  rod-elements,  defined  as  “digital  elements”,  as  shown  in 
Fig.l.  Frictionless  pins,  defined  as  “nodes”,  connect  rod-elements.  As  the  length  of  these  rod- 
elements  approaches  zero,  the  digital  chain  becomes  fully  flexible.  It  is  thus  able  to  represent  a 
one-dimensional  flexible  physical  entity  with  a  fixed  cross-section,  such  as  fibers  and  yams. 


Digital  Rod  Elements 


Frictionless  Pins 


Fig.l  Concept  of  Digital  Chain 

In  digital  element  simulation,  the  flexible  nature  of  a  digital  chain  is  conveyed  by  frictionless- 
pins,  which  connect  digital  elements.  Digital  element  size  must  be  very  small.  Otherwise, 
physicality  (in  this  case,  flexibility)  cannot  be  preserved.  The  digital  chain  is  a  digital 
representation  of  a  flexible  physical  entity.  Digital  element  length  reflects  digital  discretization 
resolution. 

The  stiffness  matrix  of  the  digital  element  is  the  same  as  the  stiffness  matrix  of  the  rod  element 
used  in  finite  element  analysis  [1-2],  which  can  be  expressed  as: 

'1  0  0-10  0 

0  A  0  0  -A  0 

EA  0  0  A  0  0  -A 

~L  -1  0  0  1  0  0 

0  -A  0  0  A  0 

0  0  -A  0  0  A 

where  E  is  the  modulus,  A  is  the  cross-section  area,  L  is  the  length  of  the  digital  element  and  A  is 
a  small  perturbation  used  to  prevent  the  singularity  of  the  global  stiffness  matrix.  A  must  be 
much  smaller  than  1 .  Generally  speaking,  it  is  on  the  order  of  1  O'6  -  1  O'12. 

2.1. 2  Contacts  between  Digital  Chains 

When  the  element  length  approaches  zero,  contact  between  two  digital  chains  can  be  represented 
by  contact  between  nodes  from  two  neighboring  chains.  See  Fig.2.  If  the  distance  between  two 
nodes  is  smaller  than  the  diameter  of  the  digital  chain,  a  contact  element  is  added  between  them. 
The  stiffness  matrix  can  be  expressed  as 


3 


Digital  chains 


Fig.2  Contact  Element  between  Two  Neighboring  Digital  Chains 


k„  0  0  -kn  0  0 

0  -Jfc,  0  0  +k,  0 

0  0  -ks  0  0  +ks 

-k„  0  0  k„  0  0 

0  +ks  0  0  -ks  0 

0  0  +ks  0  0  -ks 

where  k„  and  ks  are  the  compression  stiffness  coefficient  and  the  lateral  stiffness  coefficient, 
respectively. 

If  contact  occurs  between  two  nodes,  one  of  two  kinds  of  physical  conditions  would  exist: 
sticking  or  sliding. 

Two  digital  chains  would  stick  together  if  fjFn  >|FV| ,  where  /j.  is  the  friction  coefficient,  Fn  is  the 
compressive  force  between  two  nodes,  and  Fs  is  shear  force  between  two  digital  chains.  In  a 

rigid  contact,  displacements  of  node  i  and  node  j  are  constrained  to  be  the  same.  Thus,  the 
penalty  method  can  be  employed.  Stiffness  coefficients  k„  and  ks  are  replaced  by  a  large  penalty 
number. 

Sliding  occurs  between  two  yams  if  juFn  <|FV|.  If  this  occurs,  the  lateral  stiffness  coefficient  ks 
would  be  zero. 
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2.1.3  Yam  Assembly 

Two  digital  element  analysis  models  have  been  developed.  One  is  called  “single  chain  digital 
element  analysis”;  the  other  is  called  “multi-chain  digital  element  analysis”.  The  former  operates 
on  the  level  of  yam;  the  latter  on  the  level  of  fiber. 

Single-chain  digital  element  analysis 

In  a  single  chain  digital  element  analysis,  a  yam  is  modeled  as  a  single  digital  chain.  It  is  a 
flexible,  one-dimensional  physical  entity  with  a  constant  cross-section,  most  commonly  circular. 
Once  the  yam  is  discretized  into  digital  element  chains,  each  element  receives  an  element 
number  and  each  nock  receives  a  nodal  number.  Then,  contact  elements  are  positioned  between 
two  neighboring  digital  chains  if  they  contact  each  other.  The  element  stiffness  matrices  are 
calculated  and  a  global  stiffness  matrix  is  assembled.  Given  specific  boundary  conditions,  nodal 
displacements  and  element  stress  and  strain  can  be  derived.  Yam  paths  inside  a  fabric  are  defined 
by  nodal  positions;  yam  tensions  are  defined  by  element  stresses. 

Fig.3  shows  a  simple  twisting  process  simulated  by  single  chain  digital  element  analysis.  Fig.3- 
a  is  the  set-up  for  the  twisting  process.  Two  spring  elements  are  placed  at  the  bottom  of  the 
yams  in  order  to  maintain  yam  tension.  First,  the  top  ends  of  the  two  yams  are  moved  upward  a 
distance  of  Ah  to  achieve  initial  yam  tension.  This  provides  a  pre-twist  yam  tension.  Then,  the 
bottom  ends  of  both  spring  elements  rotate  along  a  circular  tract.  Thus,  a  twisted  yam  is 
produced.  Fig.3 -b  is  the  twisted  yam  created  by  digital  element  simulation.  One  also  can 
visualize  the  length  of  the  digital  elements  from  3-b. 


3 -a  Numerical  Model 


3-b  Twisted  Yam 


Fig.3  Digital  Model  and  Simulation  Results  for  the  Twisting  Process 


Multi-chain  digital  element  analysis 


Multi-chain  digital  element  analysis  follows  physicality  more  precisely.  In  this  model,  each  fiber 
is  modeled  as  a  digital  chain  and  each  yam  is  composed  of  many  fibers,  i.e.  a  yam  is  modeled  as 
an  assembly  of  digital  chains.  After  fibers  are  discretized  into  digital  chains  and  after  yam  is 
assembled,  the  process  used  in  the  previous  model  can  be  adopted.  Fiber  paths  are  defined  by 
nodal  positions.  Yam  is  defined  as  an  assembly  of  fibers.  Fiber  tension  and  strain  are  defined  by 
element  stress  and  strain.  Yam  tension  is  the  summation  of  fiber  tensions  within  the  yam.  Yam 
cross-section  shape  is  defined  by  fiber  arrangement  within  the  yam.  The  method  can  be  used  to 
simulate  textile  processes  and  fabric  deformation  under  both  static  and  dynamic  loads. 
Compared  to  the  single  chain  digital  element  approach,  the  multi-chain  model  can  much  more 
accurately  derive  both  yam  movement  and  cross-section  deformation.  Yam  cross-section 
deformation  influences  yarn  paths  during  forming  and  deformation  processes.  Both  yam  and 
fiber  tension  can  thus  be  predicted.  This  allows  investigation  of  fabric  failure  mechanisms  under 
external  loads. 

During  the  textile  process,  inter-fiber  compression  and  friction  play  important  roles  in  fabric 
deformation  and  failure.  Because  the  multi-chain  digital  element  approach  closely  follows 
physicality,  it  should  provide  a  much  more  accurate  result  than  the  single  chain  digital  element 
analysis  for  the  prediction  of  fabric  micro-geometry,  deformation  and  strength. 

Physically,  a  yam  is  composed  of  hundreds  or  thousands  of  fibers.  Current  computer  power 
limits  the  amount  of  chains  that  can  be  used  to  model  each  yam.  In  most  cases,  19-50  digital 
chains  are  sufficient  to  represent  yam  cross-section  geometry. 


4-a  Yam  under  Single  Side  Compression 


4-b  Deformed  Yam  4-c  Middle  Cross  Section 

Fig.4  Yam  Deformation  under  Single  Side  Compression 
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Fig.4  illustrates  the  simplest  example  of  multi-chain  digital  element  analysis.  A  yam  is  under 
single  side  compression.  It  is  modeled  as  an  assembly  of  32  digital  chains.  Fig.4-a  shows  the 
loading  condition.  Fig.4-b  is  an  isometric  view  of  the  deformed  yam.  Fig.4-c  is  the  middle 
cross-section  shape  of  the  deformed  yam. 

2.2  Numerical  Examples 

Two  numerical  examples  are  presented  in  this  section.  In  the  first  example,  a  two-dimensional 
weaving  process  is  simulated  using  both  single-  and  multi-  chain  models.  In  the  second  example, 
a  three-dimensional  braiding  process  is  simulated.  One  can  compare  the  difference  between  the 
microstructures  calculated  by  use  of  the  two  different  approaches. 

2.2.1  Two-Dimensional  Weaving  Process 

A  two-dimensional  weaving  process  was  simulated  using  the  single  chain  digital  element 
approach.  A  step-by-step  simulation  of  the  weaving  process  is  illustrated  in  Fig.5.  Tensions 
applied  to  wefts  and  warps  are  the  same  through  out  the  process.  Fig.6  is  a  woven  fabric 
generated  by  the  single  chain  digital  element  model.  Although  yams  are  compressed  from  the 
top  or  the  bottom  during  the  weaving  process,  this  model  cannot  simulate  yam  cross-section 
deformation,  as  shown  in  the  figure  6. 


Fig.5  Step-by-Step  Simulation  of  the  2-D  Weaving  Process 
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Fig.6  2-D  Woven  Fabric  Generated  by  Single  Chain  Digital  Element  Model 


Fig.7  2-D  Woven  Fabric  Generated  by  Multi-Chain  Digital  Element  Model 

Fig.7  displays  a  woven  fabric  generated  using  the  multi-chain  digital  element  model  with  the 
same  weaving  process.  Each  yam  is  an  assembly  of  19  digital  chains.  Total  cross-section  area 
of  the  19  fibers  is  equal  to  the  cross  section  area  of  the  single  yam.  Other  parameters  taken  in  the 
simulation  are  the  same.  Fig.  8  illustrates  yam  cross-section  along  different  sections  inside  the 
fabric.  Fig.8-a  is  the  cross-section  along  the  warp  direction.  Fig.8-b  is  the  five  cross-sections 
along  the  weft  direction.  The  positions  of  the  five  cross-sections  are  defined  in  Fig.8-a.  For  the 
sake  of  clarity,  cross-sections  of  fibers  that  belong  to  the  same  yam  have  the  same  color. 
Although  only  1 9  fibers  are  used  for  each  yam,  one  can  observe  that  yam  cross-section  varies 
from  one  position  to  another. 


W-5  W-4  W-3  W-2  W-1 


8-a  Side  View  along  Warp  Direction 


W-3 


W-5 


8-b  Cross-Sections  along  the  Weft  Direction 


Fig.8  Cross-Sections  of  2-D  Woven  Fabric  Generated  by 
Multi-Chain  Digital  Element  Model 
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2.2.2.  Three-Dimensional  Braiding  Process 

The  microstructure  of  the  three-dimensional  braided  preform  is  much  more  complex  than  the 
conventional  two-dimensional  textile  preform.  Yam  topology  and  geometry  inside  a  preform 
cannot  be  observed  from  the  surface  of  the  fabric.  Because  yams  are  entangled  and  because  yam 
axes  inside  most  textile  preforms  extend  as  a  three-dimensional  spatial  curve,  it  is  very  difficult 
to  investigate  detailed  yam  geometry  even  with  the  use  of  a  microscope. 

Fig.9-a  is  a  three-dimensional  braided  preform  generated  by  use  of  the  single  chain  digital 
element  model.  Fig.9-b  is  the  cross-section  of  the  preform.  Fig.lO-a  shows  a  three-dimensional 
braided  preform  generated  by  use  of  the  multi-chain  digital  element  model.  Fig.lO-b  is  the 
cross-section  of  the  preform.  Although  these  two  preforms  have  the  same  topology,  their 
respective  micro-geometries  differ.  Fig.lO-c  is  the  isometric  view  of  the  multi-chain  preform. 
Yam  cross-section  shapes  vary  from  section  to  section.  Cross-section  deformation,  in  tem, 
affects  the  spatial  paths  of  yams  and  the  volume  fraction  of  fibers.  Both  the  spatial  paths  of 
yams  and  the  volume  fraction  of  fibers  play  critical  roles  in  the  determination  of  the  mechanical 
properties  and  failure  mechanisms  of  textile  composites. 


9-a  3-D  Braided  Preform  9-b  Preform  Cross-Section 

Fig.9  3-D  Braided  Preform  Generated  by  the  Single  Chain  Digital  Element  Model 
2.3  Conclusions 

In  this  project,  digital  element  analysis  is  used  to  simulate  textile  processes.  Yet  the  value  of  the 
approach  is  much  more  extensive.  It  is  a  general  numerical  approach:  it  can  be  used  to 
investigate  other  mechanics  problems,  such  as  the  deformation  of  textile  preforms  during  the 
manufacturing  process  of  textile  composites.  Further,  it  enables  calculation  of  fiber  tension 
inside  fabrics.  As  such,  it  also  can  be  employed  to  investigate  the  strength  and  failure 
mechanisms  of  textile  fabrics  under  both  static  and  dynamic  loads.  Given  specific  boundary 
conditions  and  specific  external  loads,  the  digital  element  model  can  also  be  used  to  simulate  the 
failure  process  step  by  step  and  to  predict  the  strength  of  the  fabric. 
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10-a  3-D  Braided  Preform  10-b  Cross-Section 


10-c  Isometric  View 

Fig.  10  3-D  Braided  Preforms  Generated  by  Multi-Chain  Digital  Element  Models 


III.  HETEROGENEOUS  FINITE  ELEMENT  ANALYSIS 
3.1  Introduction 

Micro-failure  often  emanates  from  the  yam  (fiber)-matrix  interface.  It  can  be  caused  by 
interfacial  stress  c  oncentration  or  interfacial  debond.  Therefore,  it  is  important  to  predict 
interfacial  stress  concentration.  Many  attempts  have  been  undertaken  using  the  conventional 
finite  element  approach.  However,  there  have  been  two  great  obstacles: 

1 .  The  modulus  of  the  fibers  (or  yams)  can  be  50-100  times  higher  than  the  matrix.  Serious 
stress  concentration  along  the  yam-matrix  and  fiber-matrix  interfaces  therefore  exists.  It  is 
important  to  predict  this  interfacial  stress  concentration.  Use  of  the  conventional  iso- 
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parametrical  displacement  element  cannot  satisfy  the  equilibrium  conditions  along  the 
interface.  It  thus  provides  poor  prediction  performance  of  interfacial  stress  concentration. 
Significant  error  may  occur  even  with  employment  of  a  fine  mesh.  A  commonly  considered 
way  to  improve  the  convergent  rate  of  stress  is  to  employ  a  mixed  or  hybrid  finite  element 
method,  as  it  is  possible  to  construct  compatible  stress  functions.  However,  due  to  the  huge 
difference  of  moduli  between  fibers  and  matrix,  there  is  a  tangential  stress  jump  at  the 
interface  between  the  two  materials.  Even  with  the  employment  of  a  compatible  stress 
function,  the  element  still  fails  to  model  the  stress  jump,  so  it  still  may  be  unable  to  converge 
to  the  actual  stress  at  the  interface. 

2.  It  is  very  difficult  to  match  element  boundaries  to  the  yam-matrix  and  fiber  matrix  interfaces 
because  textile  composite  microstructure  is  very  complex..  The  micro-geometry  of  the  yam 
cross  section  varies  from  cross-section  to  cross-section,  so  mesh  generation  is  a  time 
consuming  and  tedious  task. 


crn"(x,y) ,  cr„'(x,y) 
T,J(x,y) 


isplacement 
={u"(x,y),  v’(x,y)} 


: Displacement 
=  {u!(x,y),  v7 (x,y)} 


Fig.  10  Interfacial  Physical  Conditions 

Fig.  10  illustrates  the  interfacial  physical  conditions  between  two  materials,  material  1  and 
material  2.  Material  1  has  a  modulus  of  Ei  and  a  Poisson  ratio  of  v/;  material  2  has  a  modulus  of 
E2  and  a  Poisson  ratio  of  v^.  The  displacement  field  in  material  1  is  defined  by  {1/ (x,y),  v'(x,y)} 
and  the  displacement  field  in  material  2  is  defined  by  {un(x,y),  vn(x,y)}.  u  denotes  displacement 
in  the  x-direction  and  v  denotes  displacement  in  the  ^-direction.  Physical  reality  at  the  interface 
includes: 

a.  continuity  of  displacement  at  the  interface, 

(Refer  to  Fig.  10.  {u(x,y),  v’(x,y)}  =  {un(x,y),  JJ(x,y)}  along  the  interface.) 

b.  discontinuity  of  normal  strains(f„),  shear  strains(y„,)  and  tangential  stress  (a,)  at  the 
interface,  and 

c.  equilibrium  conditions,  i.e.  continuity  of  normal  stress  and  shear  stress 

(<t7  =  cr'J  ,  r7,  =  t %  )  at  the  interface. 

Where  a  denotes  normal  stress,  r  denotes  shear  stress,  s  denotes  normal  strain,  y  denotes  shear 
strain  and  subscripts  n  and  t  denote  normal  direction  and  tangential  direction  at  the  interface  of 
two  materials,  respectively. 
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A  heterogeneous  element,  in  contrast  to  conventional  elements,  can  contain  one  or  more  kinds 
of  materials  with  different  moduli  and  other  mechanical  properties.  Matching  element  boundaries 
to  fiber-matrix  interfaces  is  unnecessary.  More  importantly,  heterogeneous  element  formulations 
more  realistically  reflect  interfacial  physics. 


In  this  project,  we  developed  two  heterogeneous  element  approaches.  One  is  a  displacement 
based  element  formulation;  the  other  a  mixed-form  element  formulation.  For  the  former, 
displacement  continuity  along  the  interface  is  guaranteed.  Stress  equilibrium  along  the  interface 
is  satisfied  in  a  weak  form.  Discontinuity  of  normal  strains(f„),  shear  strains(y„,)  and  tangential 
stress  (o?)  at  the  interface  can  also  be  modeled.  For  the  latter,  all  three  conditions,  continuity  of 
displacement  at  the  interface,  discontinuity  of  normal  strains,  shear  strains  and  tangential  stress 
at  die  interface  and  equilibrium  conditions,  are  strictly  satisfied. 

Numerical  results  show  the  heterogeneous  element  approach  as  much  more  accurate  than  the 
conventional  iso-parametric  and  conventional  mixed  element  approaches.  It  also  avoids 
difficulties  inherent  to  mesh-generation. 

3.2  Displacement  Based  Heterogeneous  Element 

3.2.1  Shape  Functions 


1 1  -a  6-node  heterogeneous  element  1 1  -b  Natural  Coordinates 

Fig.l  1  6-Node  Standard  Heterogeneous  Element  and  Their 
Natural  Coordinate  System  Mapping 
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Fig.l  1  shows  a  standard  6-node  heterogeneous  element,  referred  to  as  “D6”  in  this  report.  The 
element  contains  two  material  parts,  materials  /  and  II.  The  interface  is  represented  by  a  straight 
line  intersecting  two  opposite  element  sides  and  creating  the  two  side  nodes.  Each  material  part 
is  treated  as  a  sub-element  that  is  mapped  to  a  2  by  2  square  in  a  natural  coordinate  system,  the  r- 
s  system.  Shape  functions  are  denoted  as  H,K,  where  i  =  1,2,  ...,  6,  representing  element  nodal 
numbers  and  K  =  I,  II,  representing  the  material  index.  The  shape  functions  are: 


material  I 

H{  =$(l-r)(l-S) 
H[  =l(l+r)(l-*) 
Hi  =}(\  +  r)(l  +  s) 
Hi  =|(l-r)(l  +  *) 


material  II 

H"  =i(l-r)(\-s) 
H"  =i(l+r)(l-5) 

tf3"  =4(1 +'•)(!+*) 

H?  =i(l-r)(l  +  s) 


The  element  coordinates  are  interpolated  as: 


(3) 


material  I 
x= 

i= 1,2,5,6 

y=  T.H'y: 

/-I, 2, 5.6 

3.2.2  Displacement  Functions 
Degrees  of  Freedom 


material  II 


i= 6, 5,3,4 


y=  'LV'y, 


(4) 


Fig.  12  Additional  Degrees  of  Freedom 


!4 


If  iso-parametric  displacement  functions  are  applied  to  each  sub-element  shown  in  Fig.  1 1,  the 
stress  equilibrium  conditions  between  the  two  materials  cannot  be  satisfied.  In  order  to  address 
stress  equilibrium  conditions,  two  additional  degrees  of  freedom  are  introduced  to  the  sub¬ 
element  node  at  the  interface,  ce,K  and  0,K  (z  =  5,  6  and  K  =  I,  II),  where  K  is  the  material  index 

and  i  is  the  nodal  number,  a*  and  0,K  are  defined  as  derivatives  of  displacements  u  and  v 
along  the  tangential  direction  of  the  element’s  left  and  right  side  boundaries  at  the  interfacial 


node,  i.e.  af  = 


\  J 


and  0*  = 


v  dti  j 


ti  directions  are  shown  in  Fig.  12.  Since  coordinate 


mapping  is  linear  and  r  is  constant  along  the  element’s  left  and  right  sides,  a*  and  P,K  can  be 
expressed  in  terms  of  natural  coordinates,  such  as: 


K 

(duK> 

2 

hu*) 

or,  = 

1  9t  J 

II 

JM  | 

ds  J, 

rdvK> 

2 

rdvK> 

Pi  = 

{  dt  J 

i  =  L>K 

l  ds  ), 

(5) 


where  L*  is  the  side  length  of  the  element  shown  in  Fig.  12. 


Displacement  Functions 

The  element  displacements  in  material  I  and  material  II  are  interpolated  as: 


Material  1 

Material  II 

u'= 

/= 1,2,5,6  (=5,6 

/=6,5,3,4  /= 6,5 

v7  =  2>/v.,+  2>/a7 

v/7=  ZNS+ZMi'P,11 

/=1,2,5,6  <=5,6 

(=6, 5,3, 4  (=6,5 

(6) 


where  N/,  M/,  N/1,  Mi"  are  displacement  interpolation  functions  of  material  I  and  material  II. 
They  can  be  expressed  as: 


Material  I 

N‘,  =^(l  +  rXl-j)2 

Ni  =i(l  +  rX3-sXl  +  ^) 


Material  II 

Afj'-tO+'-Xs+iXi-*) 

Ml1 

<=^0  +  '-Xl-*2) 


It  can  be  proved  that  the  above  displacement  functions  satisfy  the  unit  properties,  i.e.. 
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1 .  At  any  node  j,j  =  1,2,  ...,6, 

N*  =0  if  i*j  and  N,K  =  1 
M ,K  =  0 

2.  At  interfacial  nodes  5  and  6, 


if  i  =  J 


(8-a) 


dM 


dtj 

dN[ 

dij 


—  =  0  if  i*  j  and 


m\  ...... 

— -  =  i  if  i  =  j 
dt: 


(8-b) 


-0 


In  addition,  the  continuity  of  displacement  at  the  interface  is  automatically  satisfied.  It  is  not 
difficult  to  prove  (u'(x,y),  v(x,y)}  =  fun(x,y),  v77 (x,y)}  along  the  interface. 


Stress  and  Strain 

The  strains  in  each  material  can  be  derived  by  differentiating  displacements  u  and  v  with  respect 
to  coordinates  x  and  y  as: 


s'  =<! 


xy 


8u 

dx 

dv 

dy 

du  dv 
dy  dx 


=  k  B/A5  B'j 

3x12  3x2  3x2 


u 

A7 


(9-a) 


£x 

11 

dx 

u 

II 

e  =  < 

X 

£y 

>  =  < 

dv 

dy 

du  dv 

■  =k 

3x12 

b;75 

3x2 

b!!«1 

3x2 

A"  > 

A  II 

lA6j 

Sr  j 


(9-b) 


where  U  denotes  {m,  v,  u2  v2  m3  v3  ua  v4  «5  v5  u6  v6}r,  a  vector  listing  of  nodal 

displacement  variables,  A f  denotes  vector  {arf  (if  }r  for  (/  =  5,6  and  K  =  I,  IT),  and  [b^  J  and 
kl  are  partial  strain-displacement  transformation  matrices  of  material  K  related  to  variables  in 
U  and  in  A f  respectively.  Stresses  can  then  be  derived  by  Hooke’s  Law. 


The  stresses  derived  from  displacement  interpolations  are  distributed  linearly  along  the  interface, 
and  the  strain  in  the  direction  is  linear  along  the  element’s  left  and  right  sides,  as  displayed  in 
Fig.  13.  In  order  Ic  maintain  stress  equilibrium  at  any  point  on  the  interface,  it  is  necessary  and 
sufficient  to  maintain  the  equilibrium  of  stresses  at  both  interfacial  nodes. 
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Fig.  13  Stress  and  strain 
distribution  along  interface  and 
element  boundaries 


Stress  along  the  Interface 


Fig.  14  Interfacial  co-ordinate  system 

The  interfacial  (t-n)  coordinate  system  is  shown  in  Fig.  14,  where  t  is  pointing  in  the  tangential 
direction  of  the  interface,  and  n  is  pointing  in  the  normal  interfacial  direction.  Stresses  in  the  t-n 
coordinate  system  in  each  material  K  (K  =  I,  II)  can  be  evaluated  by  using  strain-displacement 
relations,  Hooke’s  Law,  and  stress  coordinate  transformation  as  follows: 
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V 

K 

<*x 

K 

£ 

li 

s 

•< 

p±. 

II 

£y  ’ 

T'n. 

where  [T„]  is  the  stress  transformation  matrix,  [D]  is  the  elastic  matrix  and  K  is  the  material 
index  (K=I,II). 

Displacement  Compatibility  between  Heterogeneous  Elements:  Stress  Equilibrium  along 
Dissimilar  Material  Interfaces:  Element  Stiffness  Matrix 

There  are  two  stress  equilibrium  conditions  along  the  interface: 

0-,,'  =  ct"  and  tJ  =  t,"  (11) 

Since  stresses  are  linear  along  the  interface,  satisfaction  of  stress  equilibrium  conditions  at  two 
interfacial  nodes  guarantees  satisfaction  at  any  point  on  the  interface.  The  implementation  of  the 
two  conditions  described  in  eqs.  1 1  at  two  interfacial  nodes,  means  the  following  four  stress 
equations  must  be  satisfied: 

Node  5: 

[(to-)2_3][i>/][[b£,]5U  +  [b^5]5A^]=  [(TCT)2_3][l>//][[Bg]5u+  [B^5]5Af  ]  (12-a) 

2x3  3x3  3x12  3x2  2x3  3x3  3x12  3x2 

Node  6: 

[(TJwlMI  k]6U  +  k'Lu  +  kelsA?  ]  (12-b) 

2x3  3x3  3x12  3x2  2x3  3x3  3x12  3x2 

in  which  [(TjJ  denotes  the  matrix  that  contains  the  2nd  and  3rd  rows  of  matrix  [Tct  ].  Af  (K  = 
I,  II)  does  not  affect  stresses  at  node  6  and  vice  versa. 


The  above  four  equations  restrict  the  four  degrees  of  freedom;  thus,  A /  and  A  f  are  not 
independent  degrees  of  freedom  given  the  interfacial  stress  equilibrium  equations  (12-a)  and  (12- 
b)  are  considered.  Once  \\  is  determined,  A f  can  be  derived  using  Eqs.(12). 


Refer  to  Fig.  15.  In  order  to  maintain  displacement  compatibility  between  two  adjacent 
heterogeneous  elements,  e.g.,  element  A  and  element  B,  the  following  relationship  must  be 
satisfied: 


(asMaO8 

(A?r=kr 


(13) 


However,  if^)4  =  (a^)8  ,  (a^  Y  and  (a^  Y  can  be  derived  from  stress  equilibrium  conditions, 

that  are  Eqs.(12)  for  element  A  and  B,  respectively.  The  derived  (a^  Y  and  (a^  Y  can  not  equal 
each  other.  This  means  the  interfacial  stress  equilibrium  and  the  compatibility  between  two 
adjacent  heterogeneous  elements  can  not  be  satisfied  simultaneously. 
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1 5-a  Maintain  Interfacial  Equilibrium  and  Compromise  Displacement  Compatibility 


15-b  Maintain  Displacement  Compatibility  and  Compromise  Interfacial  Equilibrium 
Fig.  15  Displacement  compatibility  and  interfacial  equilibrium  conditions 


weak  form  ox  interfacial  equilibrium  conditions  is  developed  in  order  to  maintain  displacement 
continuity  between  adjacent  heterogeneous  elements.  Refer  to  Fig.  16.  For  each  element, 
instead  of  using  the  continuity  of  normal  and  shear  stresses  at  nodes  5  and  6,  an  average  of  the 
equilibrium  conditions  at  an  interfacial  node  for  the  adjacent  elements  is  used.  For  example,  on 
the  interfacial  node  between  A  and  B,  a  weak  form  of  the  equilibrium  conditions  can  be 


expressed  as: 


Through  mathematical  manipulation,  a  weak  form  of  the  equilibrium  conditions  becomes: 

Mv+Mu' +M(a;)*  +[s' KkY -[s;;E(a;T -ktM  =»  cm 

where  [sXM.[sa£.[SaI  .fsjf  and  [»!•!  E  are  derived  from 
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m; =[(T,)J_3]'(&j']k]-[D"]k]); 

k]‘=[(T.)M]'k]k,]' 

where  K  =  I,  II,  s  ~  A,  B  and  i  =  5, 6.  Assuming  (a16  J  =  (a7  f  and  (a"  f  =  (a"  Y ,  one  derives: 

+(kt +ktyiy =<>  o’) 


Fig.  16  Element  A  is  adjacent  to  elements  B  and  C 


Therefore,  (Aj  y*  and  (a"  Y  is  dependent  on  each  other.  Without  loss  of  generosity,  assume 

(a')4  is  an  independent  variable  and  (Aj7)^  is  a  dependent  variable,  which  can  be  expressed 
as: 

(A?r =(kt +(kf +kf )kf)  on 

Following  exactly  the  same  procedure,  dependent  variables  (a77  Y  is  shared  by  elements  A  and 
C,  which  can  be  expressed  as: 

(A  "  P  =  f  +  [s"  t )"  ([S„  E  u"  +  [s„  t  U  -  +  |s  +  [s^  t  )(A'  p )  (18-b) 

Dependent  variables  can  be  eliminated  in  strain-displacement  relations  and  element  stiffness 
matrices.  Substituting  Eqs.(18)  into  Eqs.(9),  the  strain-displacement  relations  are  be  expressed 
as 
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£J=\ K  0  0  B'5  B'J 


VA 

VB 

uc 

A  Ij 


(19) 


JJ 


=[kr  w  (b sr  (» ;,r  k.r] 


ir 

ui 

uc 


A  { 


6J 


where 


k  r = k  1+ tefe  x + k  X )'  [s„  x + telk  f + Is"  X  )X  x 
kr-tefef+kt^kr 

kr=kfci+kt)"[s„E'  (2o> 

ter = te  Jk  x*  k  x)'  k  x + k  f ) 
ter  -  telkf +kt )'kf ) 

The  element  stiffness  matrix  K  can  be  derived  as: 

[k]=  jiBf,  0  0  0  OB'^i 

+ j|b  ‘Y  (b rj  (b’J  k»r  k.rMkr  k)*  kt  ter  te.rl* 

II 

Stress  Boundary  Conditions 


(21) 


The  distributed  loan  applied  to  the  boundary  is  transferred  into  equivalent  nodal  forces  based 
upon  the  Virtual  Work  Principle.  Fig.  17  shows  an  example.  The  distributed  stress  in  the  x 
direction  is  applied  to  each  material  at  element  boundary  2-5-3.  The  virtual  work  done  by  the 
forces  on  virtual  displacement  Su  is: 
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8 w'  =  P  y*  q1  (r,s)8uI  (r,s)dl' 

•*2.^2 

=  — £  8ul  (r  =  l,$)cfe 

=  £  [j(l  -  sf  8u2  +{(3  - ^Xl  +  s)Su5  - |(l  -  s2  )lMc5‘8a‘  \b 


3  2  3  5  6  5 


Sw"  =  V^q"  (r,s)8u"  (r,s)dln 

II  rll 

=  —  y—  £  8u"  (r  =  1, .s)cfa 

=  ^-y-  |  [j(3  +  jX1  - 4&s  +i(l  +.*f  <5u,  +i(l -.v’  jtSa,"  /A/c,"]rf.v 

2fJL,.  V&L*.* 


2  q  Ls  ~  qU  < 
— — — Su ,  +  - — -<5m,  +- 
3  5  3  3 


(22-a) 


(22-b) 


Fig.  17  Equivalent  nodal  forces  for  the  heterogeneous  element 
based  on  the  Principle  of  Virtual  Work 

Therefore,  according  to  the  Principle  of  Virtual  Work,  the  equivalent  nodal  forces  are: 


r*I  Tl  TIJ  Tu 

(/,)s=2xi+“;  (fX 

t~,  \ » 


where  (£),  is  the  force  on  displacement  w„  and  (/«*)/  is  the  generalized  force  on  displacement 


derivative  a,  . 


The  equivalent  forces  on  dependent  variables  and  Ag  are  transferred  onto  other 
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independent  variables  based  on  the  Principle  of  Virtual  Work.  Assume  Mf  is  a  vector  listing  of 
forces  on  dependent  variables  Af .  Refer  to  Eq.(18-a).  Replacing  [sAr,[Si/f,and[s^  with 
zero  and  replacing  [sf]^,  [Sf/jf, and  [sAJ*  with  [s  A  ^ ,  [S^ ]5 , and  [sa  ^  ,  Eq.(18-a)  becomes: 

(Af )  =(kI)”'([St/]5U  +  [s/A]5(A^)  )  (24) 

The  virtual  work  by  M"  on  virtual  displacement  SA"  is: 

[m/  ]  SA1'  =  [m5"  f  ([s"  1)  '([St/f-SU  +  [st  l(s\? )  )  (25) 

Then,[M?f([s"i)  ‘is,,],  and  pi],  are  added  onto  nodal  forces 

corresponding  to  variables  U  and  A ;5  respectively. 

Displacement  Based  Heterogeneous  Element  with  Different  Interface  Locations 

The  interface  can  be  located  in  many  places  in  a  displacement  based  heterogeneous  element.  For 
the  convenience  of  mesh  generation  ,  six  elements  with  different  interface  locations  are  derived, 
as  shown  in  Fig.  18.  They  are  formulated  using  exactly  the  same  process  as  the  6-node  standard 
element;  only  the  node  mapping  sequence  from  the  element  coordinate  system  to  its  natural 
coordinate  system  is  different.  Details  can  be  found  in  reference  [2]. 

Convergence  Properties 

Displacement  functions  of  displacement  based  heterogeneous  elements  can  represent  all  rigid 
body  displacement  modes  and  represent  constant  stress  and  strain  states.  They  are  complete. 
Further,  displacement  compatibility  is  guaranteed  along  the  interface  between  adjacent 
heterogeneous  elements.  In  addition,  they  are  compatible  with  the  neighboring  iso-parametrical 
elements.  Heterogeneous  element  formulation,  therefore,  satisfies  the  condition  for  monotonic 
convergence.  Numerical  results  show  that  it  converges  at  a  much  faster  rate  than  conventional 
iso-parametric  elements. 


23 


(c)  5-node  corner  element  (d)  4-node  corner  element 


Fig.  18  2-D  displacement-based  heterogeneous  elements 
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3.3  Mixed  Heterogeneous  Elements 

As  has  been  explained,  a  displacement  based  heterogeneous  element  formulation  can  only  satisfy 
a  weak  form  of  interfacial  equilibrium  conditions.  For  this  reason,  a  second  mixed  form  of 
heterogeneous  element  is  developed.  This  formulation  enables  enforcement  of  interfacial 
equilibrium  conditions  and  compatibility  conditions. 

There  are  two  kinds  of  stress  interpolation  functions  used  in  mixed  elements: 

1.  The  u/a-c  formulation  (c  stands  for  continuity  of  stresses  between  elements),  in  which 
element  stresses  are  defined  by  nodal  stress  variables,  and  iso-parametric  interpolation 
functions; 

2.  The  u/a  formulation,  in  which  stress  variables  are  local.  Continuity  of  stress  is  not  enforced 
between  elements  but  results  from  the  solution  of  the  finite  element  if  the  mesh  is  fine 
enough. 

The  first  stress  formulation  can  represent  the  continuity  of  normal  stress  (traction  force)  and 
shear  stress  along  the  dissimilar  material  interface,  but  cannot  reflect  the  tangential  stress  jump 
along  the  interface.  The  second  approach  allows  for  the  tangential  stress  jump,  but  cannot  satisfy 
the  continuity  of  normal  stress  and  shear  stress  along  the  interface. 

The  heterogeneous  element  formulation  developed  in  this  project  modifies  the  second  approach. 
It  guarantees  the  continuity  of  normal  and  shear  stress  along  the  interface,  yet  it  also  allows  for 
the  tangential  stress  jump  between  two  materials.  In  addition,  the  displacement  function 
guarantees  compatibility  along  the  interface  and  between  elements.  Using  this  mixed  element 
approach,  the  material’s  constitutive  relations,  i.e.,  Hooke’s  Law,  is  realized  through  Reissner’s 
Variational  Principle. 

3.3.1  Mixed  Heterogeneous  Element  Formulations 

For  the  mixed  heterogeneous  formulation,  the  key  issue  is  how  to  construct  the  displacement 
and  stress  interpolation  functions.  The  stress  interpolation  should  not  be  of  too  high  an  order 
compared  to  the  displacement  interpolation,  for  the  element  could  again  behave  like  a 
displacement-based  element,  so  be  ineffective.  On  the  other  hand,  the  stress  interpolation 
should  not  be  too  low  an  order  compared  to  the  displacement  interpolation,  for  the  stress 
prediction  could  be  too  inaccurate.  Generally,  stress  interpolation  functions  are  of  a  lower 
order  than  displacement  interpolation  functions  because  actual  stresses  are  found  through 
strains  as  derivatives  of  displacements. 

Shape  Functions,  Displacement  Functions  and  Strain-displacement  relations 

Refer  to  Fig.  19.  A  standard  mixed  heterogeneous  element  contains  two  sub-elements,  each 
being  composed  from  a  material,  material  I  and  material  II ,  respectively.  Shape  functions  and 
displacement  functions  of  mixed  heterogeneous  elements  are  the  same  as  displacement  based 
heterogeneous  elements,  which  are  expressed  by  Eqs.(9)  and  can  be  written  as: 
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displacement^  w77  v77 } 
stress  =  \a" ,  any  ,  r77 } 
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Fig.  19  6-node  Standard  Mixed  Heterogeneous  Element 
with  Natural  Coordinate  System  Mapping 
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where  V  is  {w,  v,  u2  v2  u3  v3  uA  v4  u5  v5  u6  v6  a[  fi\  a'6  fil  a"  fi"  a"  ft" }  ,  a  vector 

listing  of  all  displacement  variables,  and  where  [B7]  and  [B7]  are  the  strain-displacement 
transformation  matrices  for  material  /  and  II,  respectively. 


Stress  Functions 

Stress  interpolation  functions  are  one  degree  lower  than  displacement  functions.  They  are 
defined  as: 
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where  (erx)07 ,  (crx)57 ,  and  (cr^7  denote  stresses  at  the  sub-element  center  (r=0,  5=0),  node  5 

and  node  6  for  material  I,  respectively.  Similarly,  (cr^  )0 77  ,  (crx  )577  ,  and  (cr^)/7  denote  stresses  at 
the  sub-element  center  (r=0,  5=0),  node  5  and  node  6  for  material  77,  respectively.  They  can  be 
written  as: 
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M/  and  M\l  are  the  interpolation  functions.  They  can  be  expressed  as 


Material  7 
M'  =1-5 

M's  =i(r  +  j) 
M67  =i(-r  +  5) 


Material  77 
M" =1+5 

M77  =-i(-r  +  5) 
Ml' =-i(r  +  s) 


(29) 


It  can  be  proved  that  the  above  stress  interpolation  functions  satisfy  the  unity  property  similar  to 
Eqs.(8-a). 


Eqs.(27)  can  be  expressed  in  matrix  form  as: 
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Refer  to  Fig.20.  Stress  functions  expressed  in  eqs.  (30)  are  linear  functions  along  the  interface 
and  the  boundary  of  the  element.  Based  upon  the  strain  functions  expressed  in  eqs.  (26),  both 
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normal  and  tangential  strains  are  linearly  distributed  along  the  interface.  In  addition,  tangential 
strain  distributes  linearly  along  the  left  and  right  boundaries  of  the  element. 


linear  s'1  ,<r" 


Fig.20  Stress  and  Strain  Distribution  along  the  Element  Boundary  and  Interface 
Interfacial  Stress  Equilibrium  Conditions 

Stress  along  interface  coordinates  can  be  derived  similarly.  Refer  to  Fig.  14.  Stress  along  the 
interface  can  be  written  as 
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where  [T„]  is  the  stress  transformation  matrix. 

There  are  two  equilibrium  conditions  that  need  to  be  maintained  at  the  interface  in  order  to 
reflect  physical  conditions  there: 


i  1  n 

1-  <*n  =°n  ; 

0  1  11 
2.  . 


Since  the  stress  functions  are  linear  functions  along  the  interface,  satisfaction  of  the  two 
conditions  at  two  ict erfacial  nodes  guarantees  satisfaction  at  any  point  on  the  interface. 
Implementation  of  the  two  conditions  at  both  interfacial  nodes  gives  rise  to  four  stress  equations: 
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Therefore,  the  four  stress  variables  are  dependent  variables.  They  can  be  eliminated.  Combining 
eqs.  (30),  (31)  and  (32),  stress  within  the  two  materials  can  be  expressed  in  interfacial 
coordinates  as: 
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(M5/T<TI),and(M6/Tcr,)I  consist  of  the  first  column  of  (MfT^and  (MfT^1).  (MfT^1)^ 
and  (Mf  T“*)2_3  consist  of  the  second  and  third  columns  of  (MfT'1)  and  (Mf  T^1) . 

The  formulation  of  the  mixed  heterogeneous  element  guaranteed  the  stress  equilibrium  condition 
along  the  interface  and  the  displacement  compatibility  along  the  interface  and  element  boundary . 

Distributed  Boundary  Stress  Conditions 

Using  the  virtual  work  principle,  distributed  boundary  stress  translates  to  nodal  forces.  Since 
the  mixed  heterogeneous  element  uses  the  same  displacement  functions,  the  same  procedure 
listed  in  eqs.  (22)  can  be  adopted.  Equivalent  nodal  force  resulting  from  distributed  load  can  be 
expressed  by  eqs.  (23). 

Reissner  s  Variational  Principle 

Reissner’s  Variational  Principle  enables  one  to  derive 
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where  {a}  and  {V}  are  expressed  in  Eqs.(34)  and  Eqs.(26),  {F}  is  an  equivalent  nodal  force 
vector  applied  on  displacement  {V},  and  [Q]  and  [C]  are 


[Q] = t  J  [s'  )<u + 1  J  [s"  (mu'  )<M 

1  u  (36) 

[C]  =  /  )^'\dA  + 1  ||m  J  [b"  ]dA 

1  II 

where  [S]  is  the  compliance  matrix  and  t  is  the  thickness  of  materials. 

Refer  to  Eq.  (35).  Stress  vector  {ct}  is  considered  as  internal  variable  vector.  It  is  expressed  in 
terms  of  the  element  displacement  as: 

W-Wfclm  (37) 

Eliminating  {a}  in  eq.(35),  one  can  derive 
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Mixed  Heterogeneous  Elements  with  Different  Interface  Locations 

Similar  to  the  displacement  based  heterogeneous  element,  a  total  of  six  mixed  heterogeneous 
elements  with  different  interface  locations  are  developed.  They  are  shown  in  Figs.  18. 


IV  NUMERICAL  EXAMPLES 


4.1  Numerical  Example  1 


Fig.  21  illustrates  the  first  example.  It  consists  of  two  materials.  Core  material  modulus  is  100 
times  ring  material  modulus.  Uniform  pressure  is  applied.  A  plane-stress  state  is  assumed.  Refer 
to  Fig.  21 .  Assume:  a0  =  1000;  a  =  1;  b  =  2 . 

An  analytical  solution  for  this  problem  exists.  Stress  distribution  in  material  I  is: 

<J  1  =  — <T, 
r  1 

(39) 

rre  =  0 

Stress  distribution  in  material  II  is: 
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where  cr]  is  the  radial  pressure  between  the  two  materials,  which  is: 


2b2  a, 


( l-v'  \ b 2  -  a2 )+  b2 (l  +  v“  )+  a2  (l  -  v"  ) 


-1 


(41) 


Element 

boundary 


Interface 


22-a  Coarse  Mesh 


Element 

boundary 


Interface 


22-b  Fine  Mesh 


Fig.22  Finite  Element  Mesh 


31 


Fig.23-a  Comparison  of  Radial  Stress  Derived  from 
Four  Different  Approaches  (Coarse  Mesh) 


Fig.23-a  Comparison  of  Radial  Stress  Derived  from 
Four  Different  Approaches  (Fine  Mesh) 
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The  two  finite  element  meshes  shown  in  Fig.  22  are  used  to  compare  the  accuracy  of  the 
heterogeneous  element  approach  to  conventional  iso-parametric  (homogeneous)  elements  and 
conventional  (homogeneous)  mixed  elements.  Fig.23a  compares  radial  stress  distribution  along 
the  OR  (Shown  in  Fig.22)  derived  from  the  coarse  mesh.  Fig.23b  is  derived  from  the  fine  mesh. 
Five  curves  appear  in  each  chart.  The  black  curve  is  derived  from  the  analytical  solution.  The 
red  curve  is  derived  from  iso-parametric  element  approach.  The  green  is  from  the  conventional 
mixed  element  approach.  The  blue  and  purple  curves  result  from  the  application  of  the 
heterogeneous  element  along  the  interface.  The  blue  curve  is  derived  from  an  approach  in  which 
displacement  based  heterogeneous  elements  are  employed  along  the  interface  and  iso-parametric 
elements  are  used  elsewhere.  The  purple  curve  is  derived  from  an  approach  in  which  mixed 
heterogeneous  elements  are  employed  along  the  interface  and  conventional  mixed  elements  are 
employed  elsewhere..  One  finds  that  the  blue  and  purple  curves  are  very  close  to  the  analytical 
solution.  It  shows  that  the  application  of  heterogeneous  elements  improves  accuracy 
dramatically. 

4.2  Numerical  Example  II 


I  1  I  I  Nl  I  I  I  J 


24-a  Cross-section 
♦  1  1 L*  1  0<>  I 


24-b  Unit  cell  and  Boundary  Conditions 
Fig.24  Numerical  Example  II 
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Table  1 


Material  properties  of  fibers  and  resin 


Material 

Young’s  modulus  (GPa),  E 

Poisson’s  ratio,  v 

Fiber 

325 

0.15 

Matrix 

3.45 

0.35 

Figure  24-a  shows  the  material  domain  of  the  second  numerical  example.  It  illustrates  a 
unidirectional  fibrous  reinforced  composite  cross  section.  Moduli  and  Poisson  ratios  of  fibers 
and  matrix  are  listed  in  Table  1.  Because  fibers  are  periodically  distributed  inside  the  material 
domain,  only  a  unit  cell  is  required  to  determine  stress  distribution.  The  unit  cell  geometry  and 
boundary  conditions  employed  in  the  numerical  analysis  are  shown  in  24-b,  which  is  extracted 
from  the  cross-section  of  24-a.  The  plane  strain  state  is  assumed. 


25-a  Coarse  Mesh 


25-b  Medium  Mesh 


25-c  Fine  Mesh 


Fig.25  Finite  Element  Mesh 

Three  finite  element  meshes  are  used:  coarse,  medium  and  fine.  Again,  the  equilibrium  along  the 
interface  cannot  be  enforced  using  either  the  iso-parametric  approach  or  the  conventional  mixed 
element  approach.  For  example,  26-a,  26-b  and  26-c  shows  the  interfacial  shear  stress  along  the 
interface  derived  from  conventional  iso-parametric  and  mixed  element  analysis.  One  finds  that 
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shear  stress  derived  from  the  fiber  side  differs  from  that  derived  from  the  matrix  side.  The 
difference  becomes  smaller  with  use  of  a  finer  element  mesh.  Theoretically,  the  difference  would 
approach  zero  and  the  stress  would  approach  the  actual  stress  as  the  size  of  the  element 
approaches  zero.  26-d  shows  the  results  with  the  application  of  heterogeneous  elements  along  the 
interface. 
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26-a  Isoparametric  and  mixed  Element  Approach  26-b  Isoparametric  and  mixed  Element  Approach 
(Coarse  Mesh)  (Medium  Mesh) 
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26-c  Isoparametric  and  mixed  Element  Approach  26-d  Heterogeneous  Element  Approach 
(Fine  Mesh)  (Coarse  Mesh) 

Fig.26  Interfacial  Shear  Stress  Derived  from  Various  Element  Meshes  and  Approaches 

Fig.27  compares  results  received  from  traditional  iso-parametric  element  analysis  and  results 
from  the  application  of  the  displacement  based  heterogeneous  element  along  the  interface.  There 
are  six  curves  in  Fig.27.  The  top  two  curves  are  shear  stress  distribution  along  the  interface  at 
fibers’  sides,  which  are  derived  from  the  conventional  iso-parametric  element  approach.  The  first 
top  one  shows  the  result  using  the  coarse  mesh,  the  second  the  result  using  the  fine  mesh.  The 
bottom  two  curves  are  shear  stress  distribution  along  the  interface  of  the  matrix  sides,  which  are 
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also  derived  from  the  conventional  iso-parametric  approach.  The  first  bottom  one  is  the  result 
using  the  coarse  mesh;  the  second  the  result  using  the  fine  mesh.  The  two  curves  in  the  middle 
are  results  from  the  introduction  of  heterogeneous  elements  along  the  interface.  Although  the 
coarse  mesh  is  used,  the  resulting  stress  curves  are  located  between  the  two  curves  derived  from 
the  iso-parametric  element  with  the  fine  mesh.  From  this  one  sees  the  results  from  use  of  a  coarse 
mesh  in  the  heterogeneous  element  is  more  accurate  than  use  of  a  fine  mesh  in  an  iso-parametric 
element  analysis! 

Fig.28  compares  results  received  from  conventional  mixed  element  analysis  to  results  from  the 
application  of  the  mixed  heterogeneous  element  along  the  interface.  Six  curves  are  shown.  The 
top  two  curves  show  shear  stress  distribution  along  the  interface  at  the  fibers’  sides,  which  are 
derived  from  the  conventional  mixed  element  approach.  The  first  top  one  is  the  result  using  the 
coarse  mesh;  the  second  the  result  using  the  fine  mesh.  The  bottom  two  curves  are  shear  stress 
distribution  along  the  interface  at  the  matrix  sides,  which  are  also  derived  from  the  conventional 
mixed  approach.  The  first  bottom  one  is  the  result  using  the  coarse  mesh;  the  second  the  result 
using  the  fine  mesh.  The  two  curves  in  the  middle  are  results  from  the  introduction  of 
heterogeneous  elements  along  the  interface.  Although  the  coarse  mesh  is  used,  the  resulting 
stress-curves  are  located  between  the  two  curves  derived  from  the  iso-parametric  element  with 
fine  mesh.  Again*  from this,  one  sees  the  results  from  use  of  a  coarse  mesh  in  the  heterogeneous 
element  is  more  accurate  than  use  of  a  fine  mesh  in  an  iso-parametric  element  analysis! 


Fig.27  Comparison  of  Shear  Stress  along  the  Interface  Derived  From  Displacement  Based 
Heterogeneous  Element  Approach  and  Iso-parametric  Element  Approach 
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Fig.28  Comparison  of  Shear  Stress  along  the  Interface  Derived  From  Mixed  Heterogeneous 
Element  Approach  and  Conventional  Mixed  Element  Approach 


Fig.29  Rectangular  Mesh 

As  aforementioned,  one  difficulty  with  use  of  finite  element  analysis  for  composites  is  the  time- 
consuming  mesh  generation  process.  The  heterogeneous  element  can  contain  more  than  one 
material.  There  is  no  need  to  match  the  element  boundary  to  the  interface.  Simplified  mesh 
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generation,  therefore,  can  be  used.  For  example,  the  unit  cell  used  in  the  previous  sub-section  can 
be  descretized  into  a  rectangular  element  mesh  as  shown  in  Fig.29.  The  dotted  red  curve 
represents  the  interface.  Elements  along  the  interface  are  heterogeneous.  Other  employed 
elements  are  conventional  iso-parametric.  Fig.30  shows  the  calculated  interfacial  normal  and 
shear  stress  distributions  using  both  simplified  rectangular  mesh  and  the  previous  mesh.  Again, 
results  derived  from  the  previous  mesh  and  from  the  rectangular  mesh  are  very  close. 
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30-a  Radial  Stress  along  the  Interface 


30-b  Shear  Stress  along  the  Interface 

Fig.30  Comparison  of  Interfacial  Stress  Derived  Using  Rectangular  Mesh  and  Regular  Meshes 


3.4  Summary  on  Heterogeneous  Elements 


Two  HeteK>geneous..el®rQent  formulations  are  developed:  displacement  based  and  mixed.  In  the 
former  formulation,  displacement  compatibility  is  enforced  along  the  interface  between  two 
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dissimilar  materials,  along  the  boundary  between  heterogeneous  elements,  and  along  the 
boundary  between  heterogeneous  elements  and  iso-parametric  elements.  The  equilibrium 
condition  is  satisfied  in  a  weak  form.  In  the  latter  formulation,  not  only  is  the  compatibility 
along  the  interface  of  dissimilar  materials  on  the  boundary  between  heterogeneous  elements  and 
conventional  mixed  elements  enforced,  but  the  equilibrium  condition  is  also  strictly  enforced. 

Numerical  results  show  that  replacing  a  homogeneous  element  at  the  interface  area  with  a 
heterogeneous  element  improves  accuracy  significantly.  In  addition,  the  mesh  generation  process 
is  immensely  simplified.  Because  one  heterogeneous  element  can  contain  two  materials,  there  is 
no  need  to  match  the  element  boundary  with  the  interface.  A  simple  rectangular  mesh  can  be 
adopted. 

IV.  IMPACT  OF  THE  PROJECT 

The  “Digital  element  approach”  is  a  unique  technology  for  textile  mechanics  that  has  been 
developed  through  the  research  project.  It  is  the  only  model  that  can  analyze  textile  fabric 
deformation  based  upon  fiber-scale  mechanics.  It  can  be  used  for  textile  process  design  of 
composite  structural  components,  for  deformation  simulation  of  textile  fabrics  during  the 
molding  process  and  for  the  prediction  of  textile  preform  micro-geometry,  such  as  yam  paths, 
yam  cross-section  deformation  and  fiber  paths.  Further,  it  can  be  used  for  dynamic  analysis  of 
textile  fabrics,  such  as  for  the  high  penetration  failure  mechanism  of  textile  fabrics  and  for  textile 
micro-/nano-  device  design. 

The  “Heterogeneous  element  approach”  is  a  new  element  method  developed  through  the  research 
project  specifically  for  composite  materials.  The  concept  is  that  an  element  can  consist  of  more 
than  one  material  The  formulation  of  the  heterogeneous  element  enforces  both  the  stress 
equilibrium  conditions  and  displacement  continuity  at  the  interface  between  fiber  and  matrix.  As 
a  result,  the  convergence  rate  of  the  heterogeneous  element  approach  is  much  faster  than  the 
conventional  iso-parametric  element  approach.  Two  kinds  of  heterogeneous  element 
formulations  have  been  developed:  A  displacement  based  heterogeneous  element  and  a  mixed 
heterogeneous  element.  Preliminary  tests  show  that  80-90%  less  elements  are  required  in  the 
vicinity  of  the  interface  if  the  heterogeneous  element  approach  is  employed. 


V.  PERSONNEL  SUPPORTED/ASSOCIATED: 


Principal  Investigator: 
Assistant  Research  scientist: 
Post  Doctorate  Associate: 
Ph.D.  Graduate  student: 
Master  graduate  student: 


Youqi  Wang  (Associate  Professor) 

Eric  Zhou 

Xuekun  Sun 

Chen  Tao 

Jason  Rahn 


V.  PUBLICATIONS 

1.  Guangming  Zhou,  Xuekun  Sun  and  Youqi  Wang,  “Multichain  digital  element 
analysis  in  textile  mechanics,”  Composites  Science  and  Technology,  2003. 


39 


2.  Chen  Tao,  “Heterogeneous  Element:  A  New  Finite  element  Method  for  the  Micro- 
Stress  Analysis  of  Composites,”  Ph.D.  Thesis,  Kansas  State  University,  2003. 

3.  Youqi  Wang,  Changjie  Sun  and  Xunkun  Sun,  Nicholas  J.  Pagano  "  Principles  for 
recovering  micro-stress  in  multi-level  analysis,"  ASTM  STP  1416,  Composite 
Materials:  Testing,  Design,  and  Acceptance  Criteria,  A.Zereick  and  A.T.  Nettles, 
Eds.,  Amoeban  Society  for  Testing  and  Materials  International,  2002. 

4.  Youqi  Wang,  Xuekun  Sun  and  Guangming  Zhou,  "Application  of  Digital  Element  in 
Textile  Mechanics,"  TEXCOMP-6,  Philadelphia,  2002. 

5.  Youqi  Wang  and  Chen  Tao,  "Analysis  of  fiber-matrix  interfacial  stress  analysis  by 
using  heterogeneous  finite  element,"  ICCM-12,  Paris,  July  1999. 

6.  Youqi  Wang  and  Xuekun  Sun, "  Digital  element  simulation  of  textile  processes," 
Composites  Science  and  Technology,  61  (2):2001 . 

7.  Youqi  Wang  and  Xueknu  Sun,  "  Determining  the  geometry  of  textile  preforms  using 
finite  element  analysis,"  1 5th  Technical  Conference  of  American  Society  for 
Composites,  Taxes  A&M,  Sept.  2000. 

8.  Xuekun  Sun  and  Youqi  Wang,  “Geometry  of  3-D  braided  rectangular  preform  with 
axial  yams,”  SAMPE  2001,  May,  2001. 

9.  Changjie  Sun  and  Youqi  Wang,  “Determine  the  size  of  “local  domain”  in  multi-scale 
analysis,”  ICCM  13,  Beijing,  June,  2001. 


VL  NEW  DISCOVERIES,  INVENTIONS,  OR  PATENT  DISCLOSES: 
Two  numerical  methods  are  developed: 

1 .  Digital  element  simulation  approach 

2.  Heterogeneous  finite  element  method 
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